#ifndef MLTOOL_OPERATIONS_H
#define MLTOOL_OPERATIONS_H

#include<cassert>
#include<algorithm>
#include"tags.h"
#include<type_traits>
#include"vector.h"
#include"matrix.h"
#include"lu.h"
#include"symmetric_matrix.h"
#include"triangular_matrix.h"

typedef double precision_type;

namespace math{

//任意个标量进行算术运算（+，-， ×， /）所得结果的类型
template<class ... Types>
struct result_value_type;

template<class T1>
struct result_value_type<T1>
{
	typedef T1 type;
};

template<class T1, class ... Types>
struct result_value_type<T1, Types...>
{
	typedef decltype(T1() * (typename result_value_type<Types...>::type())) type;
};

//向量×向量，向量×矩阵，矩阵×矩阵，矩阵×向量的结果类型
template<class T1, class T2>
struct category_multiplies;

template<>
struct category_multiplies<matrix_tag, matrix_tag>
{
	typedef matrix_tag type;
};

template<>
struct category_multiplies<matrix_tag, vector_tag>
{
	typedef vector_tag type;
};

template<>
struct category_multiplies<vector_tag, matrix_tag>
{
	typedef vector_tag type;
};

template<>
struct category_multiplies<vector_tag, vector_tag>
{
	typedef matrix_tag type;
};

//Category 为matrix_tag或 vector_tag
template<class Category, class T>
struct mix_type;

template<class T>
struct mix_type<matrix_tag, T>
{
	typedef matrix<T> type;
};

template<class T>
struct mix_type<vector_tag, T>
{
	typedef vector<T> type;
};

template<class T1, class T2>
struct result_type_multiplies
{
	typedef typename category_multiplies<typename T1::category, typename T2::category>::type category;
	typedef typename result_value_type<typename T1::value_type, typename T2::value_type>::type value_type;
	typedef typename mix_type<category, value_type>::type type;
};


//multiplies of two vectors.
template<class Vec1, class Vec2>
typename result_type_multiplies<Vec1, Vec2>::type
prod(const Vec1 &v1, const Vec2 &v2, vector_tag, vector_tag)
{
	typename result_type_multiplies<Vec1, Vec2>::type ans(v1.size(), v2.size());
	for(int i = 0; i < v1.size(); ++i)
		for(int j = 0; j < v2.size(); ++i)
			ans(i, j) = v1[i] * v2[j];
	return ans;
}

//multiplies of two matrixs
template<class Mat1, class Mat2>
typename result_type_multiplies<Mat1, Mat2>::type
prod(const Mat1 &m1, const Mat2 &m2, matrix_tag, matrix_tag)
{
#ifdef CHECK_DIMENSION_MATCH
	assert(m1.size2() == m2.size1());
#endif
	typename result_type_multiplies<Mat1, Mat2>::type ans(m1.size1(), m2.size2(), 0);
	for(int i = 0; i < m1.size1(); ++i)
		for(int j = 0; j < m1.size2(); ++j)
			for(int k = 0; k < m2.size2(); ++k)
				ans(i, k) += m1(i, j) * m2(j, k);
	return ans;
}

//multiplies of a matrix and a vector
template<class Mat, class Vec>
typename result_type_multiplies<Mat, Vec>::type
prod(const Mat& m, const Vec& v, matrix_tag, vector_tag)
{
#ifdef CHECK_DIMENSION_MATCH
	assert(m.size2() == v.size());
#endif
	typename result_type_multiplies<Mat, Vec>::type ans(m.size1(), 0);
	for(int i = 0; i < m.size1(); ++i)
		for(int j = 0; j < m.size2(); ++j)
			ans[i] += m(i, j) * v[j];
	return ans;
}

//multiplies of a vector and a matrix
template<class Vec, class Mat>
typename result_type_multiplies<Vec, Mat>::type
prod(const Vec&v, const Mat& m, vector_tag, matrix_tag)
{
#ifdef CHECK_DIMENSION_MATCH
	assert(v.size() == m.size1());
#endif
	typename result_type_multiplies<Vec, Mat>::type ans(m.size2(), 0);
	for(int i = 0; i < m.size1(); ++i)
		for(int j = 0; j < m.size2(); ++j)
			ans[j] += v[i] * m(i, j);
	return ans;
}

template<class T1, class T2>
typename result_type_multiplies<T1, T2>::type
prod(const T1& A, const T2& B)
{
	return prod(A, B, typename T1::category(), typename T2::category());
}

template<class T>
matrix<T> eye(size_t n)
{
	matrix<T> ans(n, n, 0);
	for(auto i = 0; i < n; ++i)
		ans(i, i) = 1;
	return ans;
}

//BLAS LEVEL 1

// y <- a * x + y
template<class T0, class Vec1, class Vec2>
inline void axpy(T0 a, const Vec1 &x, Vec2 &y)
{
#ifdef CHECK_DIMENSION_MATCH
    assert(x.size() == y.size());
#endif
    for(auto i = 0; i < x.size(); ++i)
        y[i] += a * x[i];
}

// ans = x^T * y
template<class Vec1, class Vec2>
inline auto inner_prod(const Vec1 &x, const Vec2 &y) -> decltype(x[0] * y[0])
{
#ifdef CHECK_DIMENSION_MATCH
    assert(x.size() == y.size());
#endif
    typedef typename Vec1::value_type T1;
    typedef typename Vec2::value_type T2;
    auto ans = T1() * T2();
    for(auto i = 0; i < x.size(); ++i)
        ans += x[i] * y[i];
    return ans;
}


// ans = || x ||_2
template<class Vec>
inline typename Vec::value_type norm_2(const Vec& x)
{
    typedef typename Vec::value_type T;
    auto ans = T();
    for(auto i = 0; i < x.size(); ++i)
        ans += x[i] * x[i];
    return ans;
}

// x = ax
template<class T0, class Vec>
inline void scale(Vec & x, T0 a)
{
    x *= a;
}

template<class Vec>
inline typename Vec::size_type index_norm_inf(const Vec & x)
{
#ifdef CHECK_BOUND
    assert(x.size() > 0);
#endif
    typename Vec::size_type index = 0;
    auto max_elem = abs(x[0]);
    for(auto i = 1; i < x.size(); ++i)
        if(abs(x[i]) > max_elem)
        {
            max_elem = abs(x[i]);
            index = i;
        }
    return index;
}


//BLAS LEVEL 2
//y <- alpha * Ax + beta * y
template<class Mat, class Vec1, class T, class Vec2>
inline void axpy(const Mat &A, const Vec1 &x, T beta, Vec2 &y)
{
#ifdef CHECK_DIMENSION_MATCH
    assert(A.size2() == x.size());
    assert(A.size1() == y.size());
#endif
    y *= beta;
    for(auto i = 0; i < y.size(); ++i)
        for(auto j = 0; j < A.size2(); ++j)
            y[i] += A(i,j) * x[j];
}


}

#endif
